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An analytical method is developed describing the approach to a finite-time singularity associated 
with collapse of a narrow fluid layer in an unstable Hele-Shaw flow. Under the separation of time 
scales near a bifurcation point, a long- wavelength mode entrains higher-frequency modes, as de- 
scribed by a version of Hill's equation. In the slaved dynamics, the initial-value problem is solved 
explicitly, yielding the time and analytical structure of a singularity which is associated with the 
motion of zeroes in the complex plane. This suggests a general mechanism of singularity formation 
in this system. 
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One of the most fundamental questions underlying a 
broad class of hydrodynamic phenomena is: How do 
smooth initial conditions evolve to form finite-time sin- 
gularities? Historically, two main classes of systems and 
phenomena have been investigated in this context [0. 
The first involves distributions of vorticity evolving un- 
der Euler's equation The second, dating back to clas- 
sical work of Rayleigh ||] , concerns the motion of inter- 
faces bounding masses of fluid undergoing fission [|-[n), 
or in other problems of pattern formation p2| . An is- 
sue on which considerable progress has been made is the 
relevance of self-similar solutions that may describe the 
region asymptotically close to the topology transition. 
However, there is little understanding of how such spe- 
cial solutions emerge from the large scale dynamics. 

We report here an approach to this initial-value prob- 
lem associated with a topological rearrangement of fluid 
interfaces. Under the assumptions of lubrication the- 
ory, an approximation valid for long-wavelength defor- 
mations of thin layers, the method is developed for the 
Rayleigh- Taylor instability of two-phase Hele-Shaw flow 
|p|,[13| . This lubrication approximation has been the focus 
of a considerable body of recent theoretical work on thin 
film flows and the spreading of drops Very similar, 

and even identical equations of motion describe phenom- 
ena as diverse as Marangoni convection fL4| , pattern for- 
mation in superconductors |l5f| and in biological systems 
jl6| , and oxidation of semiconductor surfaces JlTj . 

The method exploits the separation of time scales oc- 
curring close to the first instability in a system of finite 
lateral extent, where the spectrum of modes is discrete. 
As in the study of normal form expansions near convec- 
tive or lasing instabilities |l8|, the slaving of the high- 
frequency modes allows the derivation of a nonlinear evo- 
lution equation for the amplitude of the first unstable 
mode. It also allows an analytic approximation of the 
singular contribution from all other modes. 

Dynamics and Separation of Time Scales: We begin 



with the equation of motion for the half-thickness h(x, i) 
of a thin layer of fluid in Hele-Shaw flow, bounded from 
above and below by mutually immiscible fluids. The 
equation of motion in rescaled form is [^) 



h t = -d x (h [h xxx + Bh x }) 



(1) 



where the Bond number B = 2gApL 2 /a is the single 
dimensionless parameter characterizing the competition 
between the surface tension a and buoyancy associated 
with density jump Ap in a gravitational field g, and where 
L is the lateral width of the Hele-Shaw cell. An unsta- 
ble density stratification corresponds to B > The 
solution h(x, t) can also be interpreted as the height of a 
fluid layer lying at a wall, and beneath another fluid. Eq. 
|l| is also relevant to recent experiments on the pinching 
of annular rings of fluid in the Hele-Shaw cell pp| . 

As in related earlier models with B = , this PDE 
has a flux form h t + V j = arising from incompressibil- 
ity, with j = h\J the hallmark of lubrication theory. The 
characteristic velocity U ~ —VP arises from Darcy's 
law, and the pressure P is set by boundary conditions 
involving surface tension and gravity. In other contexts, 
the velocity has the more general form U ~ -/i m VP, 
such as in the spreading of drops (to = 2). Eq. |l| with 
a different physical meaning to B arises in the dynam- 
ics of the population density h of feeding herbivores [fl6[ , 
and also in the long- wavelength limit of the homogenized 
model of Type-II superconductors [^5| , with h the local 
density of vortices. 

It is sufficient for our purposes to consider Eq. ^ in a 
system of length 2tt with periodic boundary conditions. If 
h(x, t) is even and periodic, then it also describes a flow 
between two rigid walls at which the interface has 90° 
contact angle. Linearizing about a flat interface h = h 
we obtain the growth rates 



>{k) = h (Bk 2 - k 4 



(k = 0, 1, 2, 



(2) 
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The number m of unstable modes scales as V B. If B < 1 
all modes are linearly stable, whereas for 1 < B < 4 the 
mode k = 1 is unstable, while those with k > 2 remain 
damped. Moreover, if B is tuned to be slightly larger 
than unity, say B = 1 + e, then the first mode evolves 
on a time scale of order e , while the others rapidly 
equilibrate, thus being slaved to the first. 

A Contracting Flow for B = 1: Exactly at the bifur- 
cation point B = 1 there is a manifold of steady-state 
solutions, ho(x) — h(l + a cosx), for any a < 1. This 
manifold is also an attractor. Let ft, = ho + 5( x (S <C 1), 
with £ of zero mean. The linearized evolution about ho is 
Q = —ho£iC x , where Lb = d xxx + Bd x , and it preserves 
(£) =0. Consider then the norm 



T = 



with .F t = -^A: 2 (fc 2 -l)|C fc | 

k 



One may verify the inequalities T t > T t (0) exp(±Ht), 
where # = 2 min ;r {ft (x)} > @. Since T t < 0, and 
further, !F t = iff C is entirely in the null space of £i 
(and then £ is a steady state solution). One inequality 
from above then gives the bound Tt > ^t(O) e~ Ht , where 
-Ft(O) < 0. And so, JFt — * as t — ► oo, which proves 
that the nullspace of £i is attracting. Actually, ho = 
h(l + a cos for) is a steady-state for any integer fc when 
B = k 2 , but is unstable to subharmonic perturbations if 
k > 1. 

Motion along the Manifold for B > 1 : As the Bond 
number is increased slightly beyond unity, the first mode 
grows in time, but we expect that the amplitudes of 
higher modes will remain small. More generally, for 
larger B, we expect a finite number m of active modes 
including those that are linearly unstable. A natural ap- 
proach then is to partition h into low (p) and high (g) 
modes. Let V m project a periodic function onto its lower 
m modes, and write 



p + q , (V m P = V , V m q = 0) 



(3) 



Substituting this decomposition into ([!]), ignoring con- 
tributions of order q 2 , invoking slaving of higher modes, 
(dtq — 0), and integrating twice with respect to x, we 
obtain 



pq X x 



p x q x + (Pxx + Bp)q = -p t - J p + C , (4) 



where J p — pCsP is the flux associated with the lower 
modes, C is an integration constant, and / = f dx'f(x'). 
Since p is periodic in x, we find rather remarkably the 
computation of q reduced to the solution of an inhomo- 
geneous Hill's equation f2(i[| . Coupling to the partition 
constraints, V m Pt — Pt and V m q = 0, gives a complete 
set of equations to determine q and p t . 

The slaving approximation (|^) and (^) is particularly 
easy to analyze in the limit B — 1 = e — > 0, for which 



there is only one active mode, thus p = h (1 + a(t) cosx), 
and we rescale (^) with 

B = 1 + € ; t = eht ; q = eQ . (5) 

At 0(e) we obtain an inhomogeneous Ince's equation [g0| , 

(1 + a cos x) Q xx + a sin x Q x + Q 

1 2 
= a T cosx — — (1 + acosi) +C, (6) 

where C = (1 + a 2 /2)/2 by the orthogonality of p and q. 

Solvability and the Amplitude Equation: The solution 
of Eq. ^, found by variation of parameters, contains a 
secular term proportional to xsinx. Its removal is the 
solvability condition that determines a T , and yields the 
nonlinear equation of motion 



a T = | (l + V 7 !" 



(7) 



For a <C 1 this yields the exponential growth a T = a + ■ ■ ■ 
of the linear stability result (||) , but this behavior crosses 
over to a much different form in the nonlinear regime 
near pinching. Defining the function 



/(a) = o log 



(8) 



with ao = a(r = 0), we find /(&o) — f{a) — t by direct 
integration of ([?]). The singularity (or "pinch") time t p 
occurs when a / 1, so /(ao) — 1 = r p , and thus 



tp — 



/(oo) - 1 
h(B-l) 



(9) 



For ao <C 1, t p ~ log(2/ao)/h(B — 1), again consis- 
tent with the exponential growth of the linear instability. 
Near the touchdown, a(t) = 1 — (t p — t) /2 + • • • . Figure 
|l| shows excellent agreement between these asymptotic 
results and numerical studies of the lubrication PDE (|l|) 
for the pinch times t p (ao), and for the minimum height 
fe-min = h(l - a(t)). 

The correction Q is found to be 



Q{x) = A + < vl — a 2 sin x tan 1 



A- sinx 
1 — A_ cosx 



(a + cosx) In (l — 2A_ cosx + A_) 



-a — A_ cos x 

.4 2. 



(10) 



where A± are the two real zeroes of the quadratic aA 2 + 
2A + a = 0, for which A + A_ = 1, and A_ < 1. As 
a y 1, X— — >— 1, and thus within this analysis the in- 
terface curvature, through Q" '(x), develops a logarithmic 
singularity. This divergence can also be interpreted as 
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the collision on the real axis of two singularities, located 
at 7r ± iln |A_| in the complex cc-plane. 

A Spectral Cascade: The finite-time singularity can be 
understood by considering the spectrum of Q, writing 
Q = J2T=2 Q k cos k x - The recursion relation for the Qk's 
obtained from Eq. ^| has the exact solution (for k > 3) 



C 4 



k 3 -k 



A 



fc 3 - k 



(11) 



where C'± are constants determined by the inhomoge- 
neous terms. Since |A+| > 1 for a < 1, the first term in 
( |TT| ) is the secular term eliminated by the solvability con- 
dition, and thus the power-law spectrum of Q is cut off 
by an exponential factor which tends to unity as a f 1. 
Full simulations show very good agreement, to very near 
the singularity time, between the form of the correction 
function (and its spectrum) with the asymptotic result 
(|l0|). Figure ||a shows a comparison between the two 
in real space, and Fig. ^b shows the strong agreement 
between pointwise estimates — ln|A| ~ ln(Qk/Qk+i) f° r 
four wavenumbers i > 1, illustrating the collapse of the 
analyticity strip width in accord with Eq. ^. At ex- 
tremely small values of /i m i n , the slaving assumptions 
should break down, and terms such as qt cannot be ne- 
glected. Indeed, A. Bertozzi has noted that the ulti- 
mately negative divergence of Q xx , while only logarith- 
mic, is inconsistent with the existence of a single touch- 
down |pl[. Her numerical studies following /i m ; n down to 
0(1O - ^ suggest instead a saturation of the curvature. 

Bifurcation of the Singularities: Finally, we consider 
values of the Bond number well beyond the bifurca- 
tion point B = 1. Using an initial condition h = 
h (1 + acosa;), with a = 0.01, Figure 3 shows how the 
single symmetric touchdown at x — ir seen for B > 1 bi- 
furcates for B > 1.55 into two asymmetric touchdowns. 
This behavior reflects the general increase in the number 
of active modes with increasing B. While the asymp- 
totics developed for B = 1 + e are not quantitatively 
valid for e = 0(1), the slaving hypothesis with an in- 
creased number of modes leads to an appealing picture 
of how the singularities arc generated in this regime. 

The ansatz p = h (1 + a{t) cosx + b(t) cos2x) is the 
simplest allowing for two singularities, and leads to a 
spectrum of the form (lj| 
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v=l 



(fc> 1) 



(12) 



where {A„} are the zeroes of the quartic 6A 4 + a A 3 + 
2A 2 + aX + b = 0, defined by two independent quantities 
Ai,2 as Ai, A7 , A2, A^ 1 . Except when any lAi^l = 1, 
two of the four zeroes lie within the unit circle in the 
complex A plane, the other two lie outside. Elimination 
of the secular solutions associated with the latter two 
yields a t (t) and b t (t). As a and b evolve in time, the 



remaining zeroes move toward the unit circle. In the first 
quadrant of the (a, b) parameter space the line b = a — 1, 
for (1 < a < 4/3) defines those pinching configurations 
with a single touchdown (at x = 7r), while the curve 
a = [8&(1 — b)} 1 / 2 for a > 4/3 is the locus of configurations 
with two touchdowns. In the former case, Ai,2 are real (as 
are Ci^), and only one zero reaches the unit circle at the 
pinch time, thus producing only one singularity. Beyond 
point a = 4/3, b = 1/3, Ai and A2 are complex conjugates 
of each other (with C\ = = R e 1 ^) and reach the unit 
circle simultaneously, producing two singularities. 

It follows from this analysis that any Ansatz for the 
active modes contained in p generates a set of zeroes 
in the complex plane. Some of these singularities will 
move toward, although not all reach the unit circle as the 
pinch time is approached. Simulations of the full prob- 
lem agree well with this picture. The numerical studies 
indicate that the nature of the singularities depends on 
the symmetry of the touchdowns; we observe a jump dis- 
continuity in h xx for asymmetric pinching, rather than 
the divergence seen with a symmetric pinch. This is as- 
sociated with a rotation of the phase ip to tt/2. We do not 
yet understand whether this behavior may be captured 
within the slaving approximation. 

This analysis thus merges two previously separate con- 
cepts in dynamical systems described by PDEs. First is 
the coupling of slaved small scales to low-mode dynam- 
ics that recalls the reduction of a dissipative PDE to an 
inertial manifold p2[ . Second is the motion of zeroes in 
the complex plane as in the reduction of certain PDEs to 
"pole dynamics" p2] , p3[ |. 
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FIG. 3. Bifurcation diagram showing singularity locations 
versus Bond number. Insets (a) and (b) show interface evo- 
lution at B = 1.25 and B = 2.0. 



FIG. 1. Comparison between numerical solution of lubri- 
cation equation (^) and asymptotic analysis for B —* 1. (a) 
Singularity time versus initial amplitude ao from numerical 
solution of Eq. |l] (solid circles) for B = 1.05, and from 
asymptotics in Eq. ^j. (b) Minimum interface height as a 
function of time. Solid lines are the results of Eq. |^ with 
ao — 0.01,0.05,0.30 from top to bottom of figure, all with 
B — 1.05; solid circles show numerical results for those same 
initial conditions. 



FIG. 2. (a) The function Q(x) in Eq. [k] obtained from the 
asymptotic analysis (solid), compared with numerical solution 
of the full PDEs for B = 1.05 (dots). Results are for four 
times ranging from close to the initial condition to near the 
singularity time, (b) Collapse of the analyticity strip width 
as a function of time, from numerical studies. For the largest 
two values of k shown, deviations from the common curve 
arise from these amplitudes lying initially beneath machine 
precision. 
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